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Abstract 


This  report  discusses  a  technique  for  digital  filtering  by  convolution  approxi¬ 
mation  which  is  an  acceptable  compromise  between  accuracy  and  speed.  This 
technique  is  applicable  where  high  accuracy  is  not  necessary  and  where  a  digital 
computer  with  elaborate  processing  capability  is  not  available. 

Most  requirements  for  digital  filtering  can  be  developed  in  terms  of  the 
approximations  suggested  in  this  report.  By  casting  the  required  filter  in  a  form 
that  is  most  amenable  to  numerical  computation,  the  accuracy  of  approximation 
is  maximized.  If  power-of-two  accuracy  is  insufficient  there  is  a  continuous  trade¬ 
off  between  accuracy  and  speed  which  involves  range  division  and  more  bits  in 
the  approximation. 
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SECTION  I. 

Introduction 

It  is  frequently  desirable  to  be  able  to  accomplish  high-speed  digital  filtering.  The  need  arises 
for  preprocessing  of  data  for  online  calculations,  studies  of  nonstationary  signals  where  the  non- 
stationarity  is  primarily  due  to  a  certain  frequency  spectrum  and  as  an  alternative  to  detailed 
harmonic  analysis  (to  specify  three  examples)  (ref  1,  ref  3). 

A  great  deal  of  effort  has  been  expended  to  make  digital  filtering  practical  in  computation 
time,  cost,  and  accuracy  of  the  technique  of  digital  filtering.  There  are  two  general  techniques 
known  to  the  author.  These  are  the  convolution  or  Ormsby  filters  ( ref  1 )  and  the  weighted  aver¬ 
age  filters  which  include  the  recursive  filters. 

This  discussion  is  not  intended  to  review  the  entire  problem  of  digital  filtering,  but  rather  to 
explore  a  technique  for  implementing  the  convolution  technique  in  a  manner  that  is  relatively 
accurate  and  very  fast  in  computation  time  for  certain  digital  computers. 


1 


SECTION  II. 


The  Theory  of  Convolution  Filters 

The  output  of  a  linear  filter  (not  necessarily  realizable)  can  be  given  as  follows: 

F*( w)  =H(w)  F(w)  (1) 

where  F*(w)  is  the  Fourier  Transform  of  the  filter  output 
F(w)  is  the  Fourier  Transform  of  the  filter  input 
H(  w)  is  the  Fourier  Transform  of  the  filter  impulse  response 

Equation  1  has  a  time  domain  representation  as  follows: 

f*(y)  =  fh(  t)  f(t-y)  dt  (2) 

J  -oo 

where  f*(y)  is  the  filter  output 
f(t)  is  the  filter  input 
h(t)  is  the  filter  impulse  response 

A  numerical  approximation  of  equation  (2),  using  the  trapezoidal  rule  for  integration,  can  be 
written  as  follows: 

f*(yj)  =  S[h(t.)  f(ti-yj)  +h(t1+I)  f(ti+t— Tj)  ]  tl+2~tl  ]  (3) 

If  the  ti  are  equally  spaced: 

ti+i-tt  =  At  (4) 

and  therefore: 

f*(yj)  =  S[h(t,)  f(t,-yj)  +  h(t1+1)  f(t1+1-yj)  ]  (5) 

2  l 

In  general  h(ti)  is  nonzero  almost  everywhere.  However,  let  us  restrict  the  considerations  to 
the  following  case: 

h(ti)  ^  0  IiAt  —  Tmln<t|<Tmax  =  I2At  (6) 

Therefore, 

f*(yj)  =  £-2  [h(t,)  f(t,-yj)  +h(t1+1)  f(t1+1-yj)]  (7) 

2  Ii 
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This  development  is  accurate  for  an  arbitrary  f(t)  but  in  fact,  practical  considerations  limit  the 
application  to  cases  where: 


/: 


f(t)  dt  =  0 


(8) 


If  the  f(t)  does  not  fulfill  this  condition  then  write  f(t)  as  follows: 
f(t)  =  f  i  ( t )  +  A 


/: 


where  f,(t)  dt  =  0 


and 


/: 


f(t)  =  A 


Then  equation  7  can  be  rewritten  as  follows: 


(9) 


At 


I2 


f*(yj)  =  -£-2  fh(t,)  (fi(ti— yj)  +  A)  +  h(t|+i)  (fi(t|+i— yj)  +  A)  J  (10) 
2  Ii 


f*(yj)  =  a  ^-1  (h(t,)  +h(tI+1)] 

2  II 


At  L2 


~z~  2  [ h( t| )  fi(U-yi)  +  h(t,+  i)  fi(t,  +  j— yj)  ] 

2  Ii 


(11) 


or 


f*(yj)  =  fi*(yj)  +  a  r*h(t)  dt  (12) 

This  shows  that  one  can  correct  for  a  nonzero  mean  process  by  subtracting  the  mean  in  the  be¬ 
ginning  and  adding  the  mean  times  the  integral  of  the  impulse  response  of  the  filter  at  the  end. 

The  previous  development  will  suffice  for  low-pass  operations  and  the  following  for  high-pass 
operations. 

fh*(t)  =f(t)-f,*(t)  (13) 

where  fh*(t)  is  a  desired  high-pass  operation 

fi*(t)  is  the  complementary  low-pass  operation 
f(t)  is  the  unfiltered  data. 
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Therefore,  a  high-pass  operation  can  be  accomplished  by  subtracting  from  the  data  the  result  of 
passing  the  data  through  a  low-pass  filter  whose  transfer  characteristic  is  the  complement  (in 
terms  of  the  band-pass)  of  the  desired  high-pass  filter. 

Implementation  of  the  convolution  filter  in  this  form  on  a  digital  computer  involves  a  multi¬ 
plication  and  a  sum  for  each  tj  for  each  jy  The  speed  of  such  a  process  in  certain  computers 
is  often  a  severe  limitation  of  the  acceptability  of  this  filtering  technique. 
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SECTION  III. 


Approximations  to  the  Filter 

The  primary  disadvantage  of  the  previous  form  lies  in  the  time  required  for  multiplication. 
However,  it  is  obvious  that  if  the  values  of  the  filter  impulse  response  were  integral  powers  of  two, 
slow  multiply  steps  could  be  replaced  with  fast  shift  steps  in  the  digital  computer.  Therefore,  the 
question  arises;  how  well  does  the  ideal  filter  characteristic  compare  with  the  filter  that  has 
values  which  are  the  nearest  (in  some  sense)  integral  power  of  two? 

It  is  clearly  possible  to  expand  h(t)  in  a  binary  series  to  any  accuracy  desired. 

h(tm)  =  u(tm)  j  =  0, 1, 2  •  •  •  (14) 

,  _  1  )  h( tm )— 0 

where  »(t.)  -  j-  h(u)<0 

0  1  otherwise 

aj  1  j  for  each  j  corresponding  to  the 
presence  of  the  )**  power  of  two 

This  is,  of  course,  the  way  a  number  is  represented  in  a  digital  computer  and  so  implementation 
of  a  multiplication  by  expansion  in  a  binary  series  is  equivalent  to  binary  multiplication  with  N 
bits  of  accuracy.  Our  discussion  concerns  implementation  with  only  a  few  bits  of  accuracy. 

Another  point  is  that  it  is  possible  to  logically  divide  the  range  of  h(t)  in  implementing  the 
convolution  so  that  one  might  expand  the  filter  weights  as  follows: 

h(tm)  =  u(tm)  2  aj2J  h(tm)^S1  (15) 

j 

=  u(tm)  (Si  -f  2  aj2J)  S!<h(tm)^S2 

j 

=  u(tm)  (Sn  +  2  aj2*)  Sn<h(tm)^Sn  +  i 

Where  the  <{  Sj  ^  are  arbitrary  level  divisions. 

For  the  sake  of  simplicity  consider  nearest  in  the  following  sense.  Choose  nt  to  be  the  integer 
value  such  that: 

|  h  ( tj )  —  u(ti)ai2ni  |  is  a  minimum 

An  integral  square  error  criterion  has  not  been  considered.  It  may  be  that  the  approximation 
could  be  improved  thereby. 

Using  the  power-of-two  approximation  one  has: 

=  -4r  2  [u(ti)ai2n‘  fi(t, — y,)  +  u(t1+i)  a,+i  2n*+1  fi(t1  + 1 — -yj)  ]  (17) 

2  Ii 
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Using  the  logical  level  division  and  the  power-of-two  approximation  one  has: 


fi*(yj)  =  4^S[u(t.)  (Sn  +  a,2n‘)  f,(t,-y,  +  u(t,+1)  (S.  +  a1  +  12n'  +  >)  f,(t,+,-yj)] 

2  Ii  (18) 

=  £s.  2[u(t,)  fi(t,-yj)  +  u(t|  +  i)  f i ( t,  + 1 — yj ) 

2  Ii 


where 


+ 


4r  2  [u(t,)a12n«  f , ( t,-rj ) 

2  I, 


d~  u(h+i)  al  +  12n*+i  f i ( ti+ 1 — yj)  ] 


(19) 


S„  <h  ( ti )  — S„  + 1 


A  computer  program  has  been  written  for  the  IBM  7044-7094  direct  coupled  system  to  gen¬ 
erate  the  approximation  for  a  filter  of  the  following  type: 


sin  (wct) 
(wct) 


(20) 


where  wc  =  27rfc 

fc  is  the  cutoff  frequency 

The  Fourier  Transform  for  this  function  is: 

H(w)  =  1  —  w0  ^  w  ^  wc 

=  0  otherwise 

The  time  domain  graphs  for  the  function  and  its  approximation  are  illustrated  in  figure  1. 

Using  an  existing  computer  program,  SYSTRAN  (ref  2),  the  Fourier  transforms  of  the 
function  shown  in  equation  20  and  its  power-of-two  approximation  were  calculated.  The  results 
are  illustrated  in  figures  2A  and  2B. 

The  filter  accuracy  for  the  division  of  the  range  of  the  filter  values  has  not  been  evaluated. 

There  is  a  further  consideration  in  decreasing  the  time  involved  in  convolving  filters  of  the 
(sin  x)/x  type.  This  function  is  even,  therefore,  the  number  of  shift  operations  can  be  reduced 
by  half.  Further,  one  can  reduce  the  number  of  shift  operations  again  since  the  function  crosses 
some  bit  levels  more  than  once.  The  number  of  shift  operations  can  also  be  reduced  by  using 
subtract  instructions  (instead  of  add)  when  the  filter  weight  is  negative. 

n.  a,2“i  2\u(t,)  Mt.-y,)  (21) 

The  internal  sum  in  equation  21  is  for  all  values  of  fi(ti— yj)  which  are  to  be  multiplied  by  a 
particular  2Dl.  The  external  sum  is  over  all  possible  2n‘. 


0 


SECTION  IV. 


The  Accuracy  of  the  Approximate  Filter 

The  power-of-two  approximation  has  been  found  to  show  an  average  error  (relative  to  the 
exact  filter)  of  about  5%  in  the  passband  and  rolloff  range  (over  a  limited  frequency  range).  In 
the  stop-band  range  one  finds  that  the  power-of-two  filter  admits  somewhat  more  signal  energy 
than  the  exact  filter.  However,  the  approximate  filter  is  still  as  much  as  60  dB  down. 

The  rate  of  rolloff  is  primarily  determined  by  the  length  of  the  convolution  (i.e.,  Ix  and  I2). 
Using  2Y4  cycles  of  (sin  x)/x,  the  rolloff  approaches  140  dB/octave. 

One  can  generate  an  exact  expression  for  the  error  of  the  approximate  filter,  but  it  is  always 
in  the  form  of  a  series.  Therefore,  short  of  evaluating  the  series  for  an  actual  case,  the  error  can¬ 
not  be  given  a  detailed  numerical  bound. 

The  appendix  shows  the  results  of  decreasing  the  length  of  the  convolution  and  increasing 
the  interval  between  samples  in  the  convolution.  This  material  is  presented  in  the  form  of  graphs 
on  which  the  pertinent  parameters  are  listed  ( see  figures  1  through  6B ) . 
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Appendix 

Examples  of  the  Spectra  of  a  Filter  and  its  Power-of-Two  Approximation 


FIGURE  1 .  Filter  Function  and  Approximation 


9 


MAGNITUDE  (ARBITRARY  UNITS) 


FIGURE  2A.  Fourier  Spectrum  of  the  Filter 
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MAGNITUDE  (ARBITRARY  UNITS)  MAGNITUDE  (ARBITRARY  UNITS) 


FIGURE  3A.  Fourier  Spectrum  of  the  Filter 


FIGURE  3B.  Fourier  Spectrum  of  the  Filter  Approximation 
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MAGNITUDE  (ARBITRARY  UNITS)  MAGNITUDE  (ARBITRARY  UNITS) 


FIGURE  4A.  Fourier  Spectrum  of  the  Filter 


FIGURE  4B.  Fourier  Spectrum  of  the  Filter  Approximation 
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MAGNITUDE  {ARBITRARY  UNITS) 


FIGURE  5A.  Fourier  Spectrum  of  the  Filter 


FIGURE  5B.  Fourier  Spectrum  of  the  Filter  Approximation 
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MAGNITUDE  (ARBITRARY  UNITS)  MAGNITUDE  (ARBITRARY  UNITS) 


FIGURE  6A.  Fourier  Spectrum  of  the  Filter 


FIGURE  6B.  Fourier  Spectrum  of  the  Filter  Approximation 
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